rm(list=ls(all=TRUE))

path <- "" #here you need to adjust the path macro to the place on your computer
setwd(path)

library(haven)
library(dplyr)
library(ggplot2)
library(Hmisc)
library(lme4)
library(texreg)
library(effects)
library(estimatr)
library(survey)
library(sjPlot)
library(lmtest)
library(plm)
library(broom)

LevNum <- function(x){
  c(round(quantile(x, na.rm=T)[2], digits=5),
    round(mean(x, na.rm=T), digits = 5),
    round(quantile(x, na.rm=T)[4], digits=5))
}

get.or.se <- function(model) {
  broom::tidy(model) %>% 
    mutate(or = exp(estimate),
           var.diag = diag(plm::vcovHC(model, method="arellano",type="HC1",cluster="country")),
           or.se = sqrt(or^2 * var.diag),
           conf.low.or = or-1.405054*or.se, #84 percent, non overlapping ci, for 95 use 1.959947
           conf.high.or = or+1.405054*or.se) %>% 
    select(term, or, conf.low.or, conf.high.or) %>% 
    rename(estimate = or,
           conf.low = conf.low.or,
           conf.high = conf.high.or)
}

# the original data file will be freely available online (http://www.poldem.eu/). For now, please ask for the latest version from the authors.

dat <- read_dta(paste0(path, "\\chapter_wide.dta"))

######################
###### Figure 4 ######
######################

dat$iso2code <- as_factor(dat$iso2code)
dat$country <- as_factor(dat$country)

# adjusting the sample, dropping parties which have not sponsored 
# protests and cabinets under which no protest has happened.

table(dat$country)
length(unique(dat$country))


# low and high contexts

cab_cntr <-  dat %>% 
  group_by(iso2code, cabinet_name) %>%
  summarise_at(vars(aglpopw_event, wparties), funs(sum(., na.rm=TRUE)))

cab_cntr <- within(cab_cntr, prt_share <- wparties*100/aglpopw_event)
table(cab_cntr$prt_share)
cab_cntr$cabinet_name[cab_cntr$prt_share>100]
quantile(cab_cntr$prt_share)[3]
mean(cab_cntr$prt_share)

diff.cntry <- c("BE" ,"DK" ,"FI" ,"FR" ,"GR" ,"IE" ,"IT" ,"LU" ,"MT" ,"NL" ,"NO" ,"PT" ,"RO" ,"SI" ,"UK")
nodiff.cntry <- c("AT" ,"BG" ,"CY" ,"CZ" ,"EE" ,"DE" ,"HU" ,"IS" ,"LV" ,"LT" ,"PL" ,"SK" ,"ES" ,"SE" ,"CH")

cab_cntr$prt_share_cat[as.character(cab_cntr$iso2code) %in% diff.cntry] <- 0
cab_cntr$prt_share_cat[as.character(cab_cntr$iso2code) %in% nodiff.cntry] <- 1

cab_cntr$prt_share_cat <- as.factor(cab_cntr$prt_share_cat)
cab_cntr$prt_share_cat <- relevel(cab_cntr$prt_share_cat, "0")
cab_cntr$all_evts <- as.numeric(cab_cntr$aglpopw_event)

dat_upper_lev_cat <- subset(cab_cntr, select=c("iso2code",
                                               "cabinet_name",
                                               "prt_share",
                                               "prt_share_cat",
                                               "all_evts"))
dat <- merge(dat, dat_upper_lev_cat)
rm(dat_upper_lev_cat, cab_cntr)

table(dat$prt_share_cat)

# Event level regression

dat$economic <- as.factor(dat$economic)
dat$NSM <- as.factor(dat$NSM)
dat$ANSM <- as.factor(dat$ANSM)
dat$political <- as.factor(dat$political)
dat$unions <- as.factor(dat$unions)
dat$otherorg <- as.factor(dat$otherorg)
dat$petition[grepl("petition", dat$action)] <- 1
dat$petition[is.na(dat$petition)] <- 0
dat$petition <- as.factor(dat$petition)
dat$strike[grepl("strike", dat$action)] <- 1
dat$strike[is.na(dat$strike)] <- 0
dat$strike <- as.factor(dat$strike)
dat$confrontation[grepl("blockade", dat$action)] <- 1
dat$confrontation[is.na(dat$confrontation)] <- 0
dat$confrontation <- as.factor(dat$confrontation)
dat$violence[grepl("violent", dat$action)] <- 1
dat$violence[is.na(dat$violence)] <- 0
dat$violence <- as.factor(dat$violence)
dat$demo[grepl("demonstration", dat$action)] <- 1
dat$demo[is.na(dat$demo)] <- 0
dat$demo <- as.factor(dat$demo)
dat$aft_elect_prt <- as.factor(dat$aft_elect_prt)
dat$bef_elect_prt <- as.factor(dat$bef_elect_prt)
dat$country <- as.factor(dat$country)

dat$daft_parl_cycle2 <- dat$daft_parl_cycle*dat$daft_parl_cycle


dat_nodiff <- dat[dat$prt_share_cat==1,] #many parties creates a non-differentiated landscape
dat_diff <- dat[dat$prt_share_cat==0,]

m_diff_1 <- glm(parties ~ economic +
             NSM + ANSM + political + 
             aglpopw_part_all +
             unions + otherorg + petition +
             strike + confrontation + violence +
             demo + aft_elect_prt + bef_elect_prt +
               poly(all_evts, 2) +
               country,
             family = "binomial",
             weights=aglpopw_event,
             data=dat_diff)

m.eff1 <- get.or.se(m_diff_1)
m.eff1 <- m.eff1[!grepl("country", m.eff1$term),]
m.eff1$facet <- "differentiated"

m_nodiff_1 <- glm(parties ~ economic +
                  NSM + ANSM + political + 
                  aglpopw_part_all +
                  unions + otherorg + petition +
                  strike + confrontation + violence +
                  demo + aft_elect_prt + bef_elect_prt +
                    poly(all_evts, 2) +
                    country,
                family = "binomial",
                weights=aglpopw_event,
                data=dat_nodiff)

m.eff2 <- get.or.se(m_nodiff_1)
m.eff2 <- m.eff2[!grepl("country", m.eff2$term),]
m.eff2$facet <- "not differentiated"

eff.m <- rbind(m.eff1, m.eff2)
eff.m <- eff.m[!grepl("Intercept", eff.m$term),]
eff.m <- eff.m[!grepl("poly(all_evts, 2)2", eff.m$term),]
eff.m <- eff.m[!grepl("poly(all_evts, 2)1", eff.m$term),]

unique(eff.m$term)

eff.m$name[eff.m$term=="economic1"] <- "Economic"
eff.m$name[eff.m$term=="prt_share_cat1"] <- "Differentiation"
eff.m$name[eff.m$term=="NSM1"] <- "Cultural lib."
eff.m$name[eff.m$term=="ANSM1"] <- "Cultural cons."
eff.m$name[eff.m$term=="political1"] <- "Political"
eff.m$name[eff.m$term=="aglpopw_part_all"] <- "Participants"
eff.m$name[eff.m$term=="unions1"] <- "Unions"
eff.m$name[eff.m$term=="otherorg1"] <- "Other org."
eff.m$name[eff.m$term=="petition1"] <- "Petition"
eff.m$name[eff.m$term=="strike1"] <- "Strike"
eff.m$name[eff.m$term=="confrontation1"] <- "Confrontation"
eff.m$name[eff.m$term=="violence1"] <- "Violence"
eff.m$name[eff.m$term=="demo1"] <- "Demonstration"
eff.m$name[eff.m$term=="daft_parl_cycle"] <- "Electoral cycle"
eff.m$name[eff.m$term=="aft_elect_prt1"] <- "After election"
eff.m$name[eff.m$term=="bef_elect_prt1"] <- "Before election"

eff.m <- eff.m[!is.na(eff.m$name),]
eff.m$name <- factor(eff.m$name,
                     levels=c("Economic",
                              "Cultural lib.",
                              "Cultural cons.",
                              "Political",
                              "Unions",
                              "Other org.",
                              "Demonstration",
                              "Petition",
                              "Strike",
                              "Confrontation",
                              "Violence",
                              "Participants",
                              "After election",
                              "Before election"))


ggplot(data=eff.m, aes(y=estimate, x=name,
                       ymin=conf.low, ymax=conf.high,
                       group=facet,
                       color=facet)) +
  geom_pointrange(aes(shape=facet),
                  position = position_dodge2(width = 0.5, padding = 0.5),
                  show.legend=FALSE) +
  geom_hline(yintercept = 1, linetype="dotted") +
  geom_line(aes(y=conf.low-100)) + # this and the next line are for the labels
  geom_point(aes(y=conf.low-100, shape=facet), size=2.5) +
  scale_x_discrete(limits = rev(levels(eff.m$name))) +
  coord_flip(ylim=c(min(eff.m$conf.low), max(eff.m$conf.high))) +
  xlab("") + ylab("") +
  theme_bw() +
  theme(legend.title = element_blank())

ggsave(paste0(path, "f3_evt_coef.eps"),
       plot=last_plot(), dpi=300, width=7, height=5)


ggsave(paste0(path, "f3_evt_coef.png"),
       plot=last_plot(), dpi=300, width=7, height=5)


